library(Seurat)
library(tidyseurat)
library(harmony)
library(SingleR)
library(ggpubr)
library(ggrepel)
library(tidyverse)
library(data.table)
source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps')

# examine interested cytokines ------
key_cytokine <- c('Il6','Ccl20','Tslp','Flt3lg','Csf2','Tnf')

# lin22 --------
ctrl.mtx <- fread('mission/FPP/uvb/GSM5266942_C5_matrix.tsv.gz')
uvb.mtx <- fread('mission/FPP/uvb/GSM5266943_UV_matrix.tsv.gz')
vtd.mtx <- fread('mission/FPP/uvb/GSM5266944_VD_matrix.tsv.gz')

ctrl.mtx[1:5,1:5]

lin.mtx <- list(ctrl = ctrl.mtx, uvb = uvb.mtx, VtD = vtd.mtx) |>
  imap(\(x, nm)column_to_rownames(x, 'V1') |> as('sparseMatrix') |> add_name_field(nm)) |>
  reduce(RowMergeSparseMatrices)

sobj.lin <- lin.mtx |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj.lin %<>% PercentageFeatureSet('^mt-', col.name = 'mito.ratio')

sobj.lin |> VlnPlot('mito.ratio', pt.size = 0)

sobj.lin %<>% filter(mito.ratio < 5)

sobj.lin %<>% quick_process_seurat(leiden = F)

mmur <- celldex::MouseRNAseqData()

sobj.lin %<>% mark_cell_type_singler(mmur, new_label = 'mmur.main')

kc.marker <- c('Krt10','Krt14','Krt15','Krt16','Ivl','Dmkn','S100a14')

sobj.lin |> DotPlot(kc.marker, cluster.idents = T)

sobj.lin |> DimPlot(group.by = 'mmur.main')

dotty.kc <- c(15,11,18,2,0,13,6,3,5,4,1,8,7) |> sort()

sobj.lin %<>% mutate(cell.type = ifelse(seurat_clusters %in% dotty.kc,
                                      'Keratinocytes', mmur.main))

sobj.lin |> DimPlot(group.by = 'cell.type')

sobj.lin |> DimPlot(group.by = 'orig.ident')

sobj.lin |> FindAllMarkers(features = kc.marker, only.pos = T) |>
  filter(p_val_adj < .05)

sobj.lin |> write_rds('mission/FPP/uvb/lin22.rds')

sobj.lin <- read_rds('mission/FPP/uvb/lin22.rds')

## IFNL & IFNLR -------
sobj.lin |> rownames() |> str_subset('^Ifnl')

sobj.lin |>
  DotPlot('Ifnlr1', group.by = 'cell.type', split.by = 'orig.ident',
          cols = c('blue','red','green'))

ifnl.lin <- sobj.lin |>
  get_abundance_sc_long('Ifnlr1') |>
  left_join(x = sobj.lin, y = _) 

ifnl.lin |>
  ggplot(aes(orig.ident, .abundance_RNA, fill = orig.ident)) +
  geom_violin() +
  stat_summary(geom = 'crossbar', fun = 'ExpMean', width = .5) +
  facet_wrap(~cell.type)

ifnl.lin |>
  summarise(avg.expr = ExpMean(.abundance_RNA),
            pos.rate = sum(.abundance_RNA > 0) / n(),
            .by = c(orig.ident, cell.type)) |>
  ggplot(aes(orig.ident, cell.type, color = avg.expr, size = pos.rate)) +
  geom_point() +
  scale_color_distiller(palette = 'RdBu') +
  theme_pubr(legend = 'right') +
  labs(title = 'Ifnlr1 expression in mice skin treated with 300mJ/cm2 UVB',
       subtitle = 'GSE173385',
       x = 'group')

## only KC ------
kc.lin <- sobj.lin |> filter(cell.type == 'Keratinocytes')

kc.lin %<>% quick_process_seurat(skip_norm = T, leiden = F)

kc.lin |> FindMarkers()

kc.lin |> DotPlot('Trpv3', cluster.idents = T)

kc.lin |> FindAllMarkers(features = 'Trpv3', only.pos = T)

dotty.v3h <- c(1,4,7,13)

kc.lin |> mark_cell_type_singler(mmur)

kc.lin %<>% mutate(cell.type = ifelse(seurat_clusters %in% dotty.v3h,
                                      'Trpv3-hi-KC', 'Trpv3-lo-KC'),
                   subgroup = str_c(cell.type, '_', orig.ident))

kc.lin %<>% filter(seurat_clusters != 15)

kc.lin |>
  DimPlot(group.by = 'cell.type') +
  ggtitle('Keratinocyte.type')

kc.lin |>
  FindMarkers(features = c(kegg_mva, key_cytokine[key_cytokine != 'Flt3lg']), group.by = 'subgroup',
              ident.1 = 'Trpv3-hi-KC_uvb', ident.2 = 'Trpv3-hi-KC_ctrl') |>
  as_tibble(rownames = 'gene') |>
  filter(p_val_adj < .05)

kc.lin |>
  DotPlot(key_cytokine, group.by = 'subgroup')

kc.lin |>
  DotPlot(kegg_mva, group.by = 'subgroup') +
  RotatedAxis()

kc.lin |> VlnPlot(late.kc, group.by = 'cell.type', pt.size = 0)

kc.lin |> DotPlot(c(late.kc, 'Trpv3'), group.by = 'seurat_clusters',
                  cluster.idents = T, cols = 'RdYlBu')

# jain24 -------
ja24.path <- list.files('mission/FPP/uvb/jain24/', '\\d', full.names = T)

ja24.mex <- tibble(path = ja24.path) |>
  mutate(type = str_extract(path, 'barc|feat|matr'),
         sample = str_extract(path, '\\d(?=.gz)'),
         group = case_match(sample, c('1','5') ~ '0hr',
                      c('2','6') ~ '4hr',
                      c('3','7') ~ '24hr',
                      .default = '48hr'),
         sample = str_c(group, '.', sample)) |>
  pivot_wider(names_from = type, values_from = path) |>
  pmap(\(sample, group, barc, feat, matr)ReadMtx(matr, barc, feat) |>
         add_name_field(sample), .progress = T)

ja24.mtx <- RowMergeSparseMatrices(ja24.mex[[1]], ja24.mex[-1])

sobj.j24 <- ja24.mtx |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj.j24 %<>% PercentageFeatureSet('^mt-', col.name = 'mito.ratio')

sobj.j24 |> VlnPlot('mito.ratio', pt.size = 0)

sobj.j24 %<>% mutate(group = str_extract(orig.ident, '.+hr') |>
                       fct_relevel('0hr','4hr'))

sobj.j24 |> VlnPlot('mito.ratio', pt.size = 0, group.by = 'group')

sobj.j24 %<>% filter(mito.ratio < 5)

sobj.j24 %<>% quick_process_seurat(leiden = F)

sobj.j24 %<>% mark_cell_type_singler(mmur, new_label = 'mmur.main')

sobj.j24 |> DotPlot(kc.marker)

sobj.j24 |> DimPlot(group.by = 'mmur.main', cols = 'Paired')

sobj.j24 %<>% mutate(cell.type = ifelse(seurat_clusters == 21,
                                        'Keratinocytes', mmur.main))

sobj.j24 |> DimPlot(group.by = 'cell.type', cols = 'Paired')

sobj.j24 |> write_rds('mission/FPP/uvb/jain24.rds')

## ~200 KC -------
kc.j24 <- sobj.j24 |> filter(seurat_clusters == 21)

kc.j24 |> ggplot(aes(group)) + geom_bar()

kc.j24 %<>% quick_process_seurat(leiden = F, skip_norm = T)

kc.j24 |> DotPlot('Trpv3')

kc.j24 |> DimPlot(group.by = 'group')

kc.j24 %<>% mutate(KC.type = ifelse(seurat_clusters == 3,
                                      'Trpv3-hi-KC','Trpv3-lo-KC'),
                   subgroup = str_c(KC.type, '_', group))

kc.j24 |> DimPlot(group.by = 'KC.type')

kc.j24 |>
  DotPlot(key_cytokine, group.by = 'group') +
  ggtitle('Total KC')

kc.j24 |> 
  filter(str_detect(KC.type, 'hi')) |>
  DotPlot(key_cytokine, group.by = 'group') +
  ggtitle('Trpv3-hi-KC')

kc.j24 |>
  DotPlot(kegg_mva, group.by = 'group') +
  ggtitle('Total KC') +
  RotatedAxis()

kc.j24 |> 
  filter(str_detect(KC.type, 'hi')) |>
  DotPlot(kegg_mva, group.by = 'group') +
  ggtitle('Trpv3-hi-KC') +
  RotatedAxis()

late.kc <- c('Krt1','Krt10','Lor','Ivl','Tgm1','Flg')

kc.j24 |> DotPlot(late.kc)

